{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "M = sy.MatrixSymbol('M', 3, 3)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}a & 2 & 2 a\\\\2 & 1 & b\\\\b & a & 3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[a, 2, 2*a],\n",
       "[2, 1,   b],\n",
       "[b, a,   3]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import Symbol\n",
    "a = Symbol('a')\n",
    "b = Symbol('b')\n",
    "A = sy.Matrix(3,3,[a,2,b,2,1,a,2*a,b,3])\n",
    "B=A.transpose()\n",
    "\n",
    "B[1,1]\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle a + 4$"
      ],
      "text/plain": [
       "a + 4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sy.trace(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'unit_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_21756/133189190.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     40\u001b[0m \u001b[0mT_B_C\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtransformation_matrix_qt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'q_B_C_x'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_y'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_z'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'q_B_C_w'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m't_B_C_x'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m't_B_C_y'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m't_B_C_z'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     41\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 42\u001b[0;31m \u001b[0mA\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mT_W_B\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mT_B_C\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0munit_x\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     44\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'unit_x' is not defined"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def transformation_matrix_qt(q,t):\n",
    "    qx=Symbol(q[0])\n",
    "    qy=Symbol(q[1])\n",
    "    qz=Symbol(q[2])\n",
    "    qw=Symbol(q[3])\n",
    "\n",
    "    px=Symbol(t[0])\n",
    "    py=Symbol(t[1])\n",
    "    pz=Symbol(t[2])\n",
    "\n",
    "    \n",
    "    \n",
    "    qx_2=Symbol(q[0]+\"*\"+q[0])\n",
    "    qy_2=Symbol(q[1]+\"*\"+q[1])\n",
    "    qz_2=Symbol(q[2]+\"*\"+q[2])\n",
    "    \n",
    "    T=sy.Matrix(4,4,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])\n",
    "     \n",
    "    T[0,0]=1-2*qy_2-2*qz_2\n",
    "    T[1,1]=1-2*qx_2-2*qz_2\n",
    "    T[2,2]=1-2*qx_2-2*qy_2\n",
    "\n",
    "    T[0,1]=2*(qx*qy+qw*qz)\n",
    "    T[2,0]=2*(qy*qz+qw*qx)\n",
    "    T[1,2]=2*(qx*qz+qw*qy)\n",
    "\n",
    "    T[0,2]=2*(qx*qz-qw*qy)\n",
    "    T[1,0]=2*(qx*qy-qw*qz)\n",
    "    T[2,1]=2*(qy*qz-qw*qx)\n",
    "\n",
    "    T[0,3]=px\n",
    "    T[1,3]=py\n",
    "    T[2,3]=pz\n",
    "    T[3,3]=1\n",
    "    \n",
    "    return T\n",
    "\n",
    "T_W_B=transformation_matrix_qt(['q_x','q_y','q_z','q_w'],['p_x','p_y','p_z'])\n",
    "T_B_C=transformation_matrix_qt(['q_B_C_x','q_B_C_y','q_B_C_z','q_B_C_w'],['t_B_C_x','t_B_C_y','t_B_C_z'])\n",
    "\n",
    "A=T_W_B*T_B_C*unit_x\n",
    "\n",
    "for i in range(A.shape[0]):\n",
    "    for j in range(A.shape[1]):\n",
    "        print(\"P[\"+str(i)+\",\"+str(j)+\"]=\\n\"+str(A[i,j]))\n",
    "        \n",
    "        print(\"   \")\n",
    "        print(\"  \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$\\left[\\begin{matrix}\\frac{\\left(p_{T x} - p_{x}\\right) \\cos{\\left(\\theta \\right)} + \\left(p_{T y} - p_{y}\\right) \\sin{\\left(\\theta \\right)}}{\\sqrt{\\left|{p_{T x} - p_{x}}\\right|^{2} + \\left|{p_{T y} - p_{y}}\\right|^{2} + \\left|{p_{T z} - p_{z}}\\right|^{2}}}\\end{matrix}\\right]$\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "\n",
    "def transformation_matrix_dh(theta,alpha,a_,d_):\n",
    "    st=Symbol(\"sin(\"+theta+\")\")\n",
    "    ct=Symbol(\"cos(\"+theta+\")\")\n",
    "    sa=Symbol(\"sin(\"+alpha+\")\")\n",
    "    ca=Symbol(\"cos(\"+alpha+\")\")\n",
    "    \n",
    "    st=sy.sin(theta)\n",
    "    ct=sy.cos(theta)\n",
    "    sa=sy.sin(alpha)\n",
    "    ca=sy.cos(alpha)\n",
    "    \n",
    "    a=Symbol(a_)\n",
    "    d=Symbol(d_)\n",
    "    \n",
    "    T=sy.Matrix(4,4,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])\n",
    "    \n",
    "    T[0,0]=ct\n",
    "    T[1,1]=ct*ca\n",
    "    T[2,2]=ca\n",
    "\n",
    "    T[0,1]=-st*ca\n",
    "    T[2,0]=0\n",
    "    T[1,2]=-ct*sa\n",
    "\n",
    "    T[0,2]=st*sa\n",
    "    T[1,0]=st\n",
    "    T[2,1]=sa\n",
    "\n",
    "    T[0,3]=a*ct\n",
    "    T[1,3]=a*st\n",
    "    T[2,3]=d\n",
    "    T[3,3]=1\n",
    "    \n",
    "    return T\n",
    "\n",
    "trans=transformation_matrix_dh(\"theta\",\"alpha\",\"a\",\"d\")\n",
    "    \n",
    "unit_x=sy.Matrix(4,1,[1,0,0,1])\n",
    "unit_y=sy.Matrix(4,1,[0,1,0,1])\n",
    "unit_z=sy.Matrix(4,1,[0,0,1,1])\n",
    "\n",
    "\n",
    "\n",
    "trans.subs('alpha',0)\n",
    "\n",
    "unitX=sy.Matrix(3,1,[1,0,0])\n",
    "\n",
    "px,py,pz=sy.symbols(['p_x','p_y','p_z'])\n",
    "ptx,pty,ptz=sy.symbols(['p_T_x','p_T_y','p_T_z'])\n",
    "pb=sy.Matrix(3,1,[px,py,pz])\n",
    "pt=sy.Matrix(3,1,[ptx,pty,ptz])\n",
    "\n",
    "vtb=pt-pb\n",
    "\n",
    "dist=vtb.norm()\n",
    "\n",
    "vc=trans[:3,:3]*(unitX)\n",
    "\n",
    "\n",
    "print(\"$\"+sy.latex(vc.transpose()*vtb/dist)+\"$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\left[\\begin{matrix}\\frac{\\left(p_{T x} - p_{x}\\right) \\cos{\\left (\\theta \\right )} + \\left(p_{T y} - p_{y}\\right) \\sin{\\left (\\theta \\right )}}{\\sqrt{\\left|{p_{T x} - p_{x}}\\right|^{2} + \\left|{p_{T y} - p_{y}}\\right|^{2} + \\left|{p_{T z} - p_{z}}\\right|^{2}}}\\end{matrix}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7febcb040790>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGiCAYAAADgPBIcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwjElEQVR4nO3de3RU5b3/8c8QJkMCSYAEMhMJISKc1sbqEtQQa0O0gVCFqi0V86smRwQVolKgVuqFBLmoB6i/hQcvRxq1Hg7Y31EXVCsECygH0JwUlwGOXCp3EigBEiAwGZLn90dORsYEDJLJJHner7Vmmf3sZ/Z8v+wRPtmzZ2+HMcYIAADAMp1CXQAAAEAoEIIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJWCGoI+/vhjjRw5UgkJCXI4HHrvvfcC1htjlJ+fr4SEBEVERGjo0KHasmVLwByv16uHH35YcXFx6tq1q0aNGqX9+/cHs2wAAGCBoIagU6dO6eqrr9aLL77Y5Prnn39e8+fP14svvqji4mK53W5lZmbqxIkT/jmTJk3Su+++qyVLlmjdunU6efKkbrvtNtXW1gazdAAA0ME5WusGqg6HQ++++65uv/12SfVHgRISEjRp0iT99re/lVR/1Cc+Pl7PPfecHnjgAVVWVqpXr1764x//qLvuukuSdPDgQSUmJuqDDz7Q8OHDW6N0AADQAXUO1Qvv2rVL5eXlGjZsmH/M5XIpPT1d69ev1wMPPKCSkhL5fL6AOQkJCUpJSdH69evPG4K8Xq+8Xq9/ua6uTkePHlVsbKwcDkfwmgIAAC3GGKMTJ04oISFBnTq1/IdXIQtB5eXlkqT4+PiA8fj4eO3Zs8c/Jzw8XD169Gg0p+H5TZkzZ44KCgpauGIAABAK+/btU58+fVp8uyELQQ2+eWTGGPOtR2u+bc60adM0efJk/3JlZaX69u2r7du3q2fPnpdWcDvi8/m0evVqZWRkyOl0hrqcVkPf9G0D+qZvGxw9elQDBw5UVFRUULYfshDkdrsl1R/t8Xg8/vHDhw/7jw653W7V1NTo2LFjAUeDDh8+rLS0tPNu2+VyyeVyNRrv2bOnYmNjW6qFNs/n8ykyMlKxsbFW/U9D3/RtA/qmb5sE61SWkF0nKDk5WW63W0VFRf6xmpoarV271h9wBg0aJKfTGTCnrKxMmzdvvmAIAgAA+DZBPRJ08uRJ7dy507+8a9cuff755+rZs6f69u2rSZMmafbs2RowYIAGDBig2bNnKzIyUtnZ2ZKkmJgYjR07VlOmTFFsbKx69uypqVOn6qqrrtJPfvKTYJYOAAA6uKCGoP/+7/9WRkaGf7nhPJ2cnBy9/vrreuyxx3T69GlNmDBBx44d0w033KCVK1cGfPb3+9//Xp07d9Yvf/lLnT59Wrfccotef/11hYWFBbN0AADQwQU1BA0dOlQXugyRw+FQfn6+8vPzzzunS5cuWrBggRYsWBCECgEAgK24dxgAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArBTyENSvXz85HI5Gj4kTJ0qScnNzG61LTU0NcdUAAKC96xzqAoqLi1VbW+tf3rx5szIzMzV69Gj/WFZWlgoLC/3L4eHhrVojAADoeEIegnr16hWw/Oyzz6p///5KT0/3j7lcLrnd7tYuDQAAdGAhD0Hnqqmp0VtvvaXJkyfL4XD4x9esWaPevXure/fuSk9P16xZs9S7d+/zbsfr9crr9fqXq6qqJEk+n08+ny94DbQxDb3a1LNE3/RtB/qmbxsEu1+HMcYE9RUuwttvv63s7Gzt3btXCQkJkqSlS5eqW7duSkpK0q5du/TUU0/p7NmzKikpkcvlanI7+fn5KigoaDS+ePFiRUZGBrUHAADQMqqrq5Wdna3KykpFR0e3+PbbVAgaPny4wsPDtXz58vPOKSsrU1JSkpYsWaI777yzyTlNHQlKTExUWVmZYmNjW7zutsrn86moqEiZmZlyOp2hLqfV0Dd924C+6dsGFRUV8ng8QQtBbebjsD179mjVqlV65513LjjP4/EoKSlJO3bsOO8cl8vV5FEip9Np1ZunAX3bhb7tQt92sa3vYPca8q/INygsLFTv3r116623XnBeRUWF9u3bJ4/H00qVAQCAjqhNhKC6ujoVFhYqJydHnTt/fXDq5MmTmjp1qjZs2KDdu3drzZo1GjlypOLi4nTHHXeEsGIAANDetYmPw1atWqW9e/fqvvvuCxgPCwtTaWmp3nzzTR0/flwej0cZGRlaunSpoqKiQlQtAADoCNpECBo2bJiaOj87IiJCK1asCEFFAACgo2sTH4cBAAC0NkIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFgp5CEoPz9fDocj4OF2u/3rjTHKz89XQkKCIiIiNHToUG3ZsiWEFQMAgI4g5CFIkn7wgx+orKzM/ygtLfWve/755zV//ny9+OKLKi4ultvtVmZmpk6cOBHCigEAQHvXJkJQ586d5Xa7/Y9evXpJqj8K9MILL+iJJ57QnXfeqZSUFL3xxhuqrq7W4sWLQ1w1AABozzqHugBJ2rFjhxISEuRyuXTDDTdo9uzZuvzyy7Vr1y6Vl5dr2LBh/rkul0vp6elav369HnjggSa35/V65fV6/ctVVVWSJJ/PJ5/PF9xm2pCGXm3qWaJv+rYDfdO3DYLdr8MYY4L6Ct/iL3/5i6qrqzVw4EAdOnRIM2fO1JdffqktW7Zo27ZtuvHGG3XgwAElJCT4nzN+/Hjt2bNHK1asaHKb+fn5KigoaDS+ePFiRUZGBq0XAADQcqqrq5Wdna3KykpFR0e3+PZDHoK+6dSpU+rfv78ee+wxpaam6sYbb9TBgwfl8Xj8c8aNG6d9+/bpww8/bHIbTR0JSkxMVFlZmWJjY4PeQ1vh8/lUVFSkzMxMOZ3OUJfTauibvm1A3/Rtg4qKCnk8nqCFoDbxcdi5unbtqquuuko7duzQ7bffLkkqLy8PCEGHDx9WfHz8ebfhcrnkcrkajTudTqvePA3o2y70bRf6tottfQe71zZxYvS5vF6v/ud//kcej0fJyclyu90qKiryr6+pqdHatWuVlpYWwioBAEB7F/IjQVOnTtXIkSPVt29fHT58WDNnzlRVVZVycnLkcDg0adIkzZ49WwMGDNCAAQM0e/ZsRUZGKjs7O9SlAwCAdizkIWj//v26++67deTIEfXq1UupqanauHGjkpKSJEmPPfaYTp8+rQkTJujYsWO64YYbtHLlSkVFRYW4cgAA0J6FPAQtWbLkgusdDofy8/OVn5/fOgUBAAArtLlzggAAAFoDIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAVgp5CJozZ46uu+46RUVFqXfv3rr99tu1bdu2gDm5ublyOBwBj9TU1BBVDAAAOoKQh6C1a9dq4sSJ2rhxo4qKinT27FkNGzZMp06dCpiXlZWlsrIy/+ODDz4IUcUAAKAj6BzqAj788MOA5cLCQvXu3VslJSX68Y9/7B93uVxyu93N2qbX65XX6/UvV1VVSZJ8Pp98Pl8LVN0+NPRqU88SfdO3Heibvm0Q7H4dxhgT1Fe4SDt37tSAAQNUWlqqlJQUSfUfh7333nsKDw9X9+7dlZ6erlmzZql3795NbiM/P18FBQWNxhcvXqzIyMig1g8AAFpGdXW1srOzVVlZqejo6BbffpsKQcYY/exnP9OxY8f0ySef+MeXLl2qbt26KSkpSbt27dJTTz2ls2fPqqSkRC6Xq9F2mjoSlJiYqLKyMsXGxrZKL22Bz+dTUVGRMjMz5XQ6Q11Oq6Fv+rYBfdO3DSoqKuTxeIIWgkL+cdi58vLy9MUXX2jdunUB43fddZf/55SUFA0ePFhJSUl6//33deeddzbajsvlajIcOZ1Oq948DejbLvRtF/q2i219B7vXNhOCHn74YS1btkwff/yx+vTpc8G5Ho9HSUlJ2rFjRytVBwAAOpqQhyBjjB5++GG9++67WrNmjZKTk7/1ORUVFdq3b588Hk8rVAgAADqikH9FfuLEiXrrrbe0ePFiRUVFqby8XOXl5Tp9+rQk6eTJk5o6dao2bNig3bt3a82aNRo5cqTi4uJ0xx13hLh6AADQXoX8SNBLL70kSRo6dGjAeGFhoXJzcxUWFqbS0lK9+eabOn78uDwejzIyMrR06VJFRUWFoGIAANARhDwEfduX0yIiIrRixYpWqgYAANgi5B+HAQAAhAIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArNQ51AUAQMicOCH9139Jp09LYWFSTIyUlFT/cDhCXR2AICMEAei4jJH27ZM+/VTatEn64gvpF7+QcnPr1+/bJ40Y0fh58fHSr34lPfmk1L17a1YMoBURggB0LFVV0ttvS6tXS2vXSgcOBK5PTPw6BPXoIV1zjRQRIdXVSUePSrt3S4cOSfPm1YenFSukyMhWbgJAayAEAWjfamvlOn786+UzZ6Rx475eDguTfvhDafDg+sBz441fr/N46o8QncvrlVaulEaNktatk/70JyknJ5gdAAgRQhCA9scYaf16ackSdf5//0+DevWSsrPr1/XuLd13n3TZZdLQoVJq6sUdyXG5pJEjpaws6cMPpWPHgtICgNAjBAFoP/bskQoLpT/+UfrqK0mSQ1J0dXX9ESCns37eokWX9jrGSD/4gbRmjXTLLZe2LQBtFiEIQPvw619L//f/1gcUSerWTXV33KkvbxqmnT1cii2rVmrXbgrr1ALf6nI4pLlzpd/9TurZ89K3B6BN4jpBANqmU6fqj+40SEmpD0A33yy99ZaKPvpcN37/Xt2xO0bG6dR9bxTrR8/9VR9uLvvur7lmTf3RpgYEIKBDazchaOHChUpOTlaXLl00aNAgffLJJ6EuCUAwVFZKzzxTf62eP/zh6/HsbGnHDumjj/Th1Tdr/DtfqqzyTMBTyyvP6KG3/nbxQejsWWnmzPqPvkaM4DwgwBLtIgQtXbpUkyZN0hNPPKFNmzbppptu0ogRI7R3795QlwagpVRXS88+K/XrJz39tFRRIb333tfrIyKkK65QbZ1RwfKtMk1somGsYPlW1dY1NaMJpaXSj34kPfVU/dfkr79e6tLl0noB0C60i3OC5s+fr7Fjx+r++++XJL3wwgtasWKFXnrpJc2ZM6fRfK/XK6/X61+uqqqSJPl8Pvl8vtYpug1o6NWmniX6bnd9GyPH4sUKe+IJOQ4erB/63vdUO22azC9/KX2jn892HdXRk6flCqtfdnUyAf+VpKMnT2vjzsO6PvkCH2cdP65Oc+ao04IFcpw9KxMTo9r582V+9av6c4La+J9ju93fl4i+7ew7WBzGmGb+uhQaNTU1ioyM1J/+9Cfdcccd/vFHH31Un3/+udauXdvoOfn5+SooKGg0vnjxYkVy0TOgTfnhyy8r+cMPJUnVvXrpf/7P/9H+m26qv75PkHT5xz80dPJkuU6ckCQdTE1V6f3360xcXNBeE8DFq66uVnZ2tiorKxUdHd3i22/zR4KOHDmi2tpaxcfHB4zHx8ervLy8yedMmzZNkydP9i9XVVUpMTFRGRkZio2NDWq9bYnP51NRUZEyMzPlbPjqsAXou5317fHIbNiguqlT5Xz0Uf2wSxf98ALTP9t1VPe9UexfdnUyemZwnZ76707y1n39zbA/5FwXeCSork7q9L9nABijsD/8Qeb4cdU+/7x6DR+um1u4rWBrt/v7EtG3XX1XVFQEdfttPgQ1cHzjZobGmEZjDVwul1wuV6Nxp9Np1ZunAX3bpc33fexY/e0osrLql6+/Xtq3T2FRUWrOsZ/UK3qrZ7cIlVeeCTgvyFvnkLfWIYckd0wXpV7Ru/7r8keOSC+9JL35Zv3rNnzj609/kuLj1blzu/lrsEltfn8HCX3bIdi9tvkTo+Pi4hQWFtboqM/hw4cbHR0C0Mbt3y+lpUk/+5m0YcPX41FRzd5EWCeHpo+8UlL9hRLP1bA8feSVCvtbSf2VoxMT60+03rkz8Ntml10mtfMABODStPkQFB4erkGDBqmoqChgvKioSGlpaSGqCsBF+8c/6m9j8eWX9Xdpj4j4zpvKSvHopV9dK3dM4Le4krsYve/YpKx7fipdd1391aXPnJEGDZIWL5YeffQSmwDQkbSLX4MmT56se+65R4MHD9aQIUP06quvau/evXrwwQdDXRqA5jBGuvde6e9/l5KT6y9K2LfvJW0yK8WjzCvd2ri9XEe2faY/5Fyn1Finwjx31n+zKzxcGj1aeuih+qNP5/n4HIC92kUIuuuuu1RRUaEZM2aorKxMKSkp+uCDD5SUlBTq0gA0x8cf19+M1OWS/vznSw5A2r9fev99hS1friHHj2v5b36j65N7KszprL+9Rp8+9RdXtOiLEAAuXrsIQZI0YcIETZgwIdRlAPgu1qyp/+/Pfy5deeXFP//wYem//qt+O6tWSVu3+ld1ktQlN/fruc89dwmFArBJuwlBANqxhuvvbNtW/9HY+T6aOnNG2r1b2r5d+slPpIbrev3ud4F3hnc4pNRU6bbb5PvpT3Xm3Pt9AUAzEYIABN+999Z/TX3evK8D0LJl9eGmU6f6W2YcPRp4z66NG6Ubbqj/+aab6r/eftNN9ff3ysj4+qvuPl/gTU8BoJkIQQCCLypK2rRJOveaH5s3S1u2NJ7brZvUv790+vTXY/feK+XkBL9OAFYhBAFoHd+86Fl2dv1HWsbUnzDds2f9V+d79mz8cRnf7AIQBIQgAKHRr1/9AwBCpM1fLBEAACAYOBIE2KiuVtqzXjp5SOoWLyWlSZ2Cd9d2AGiLCEGAbbYukz78rVR18Oux6AQp6znpylGhqwsAWhkfhwE22bpMevvewAAkSVVl9eNbl4WmLgAIAUIQYIu62vojQDJNrPzfsQ8fr58HABYgBAG22LO+8RGgAEaqOlA/DwAsQAgCbHHyUMvOA4B2jhAE2KJbfMvOA4B2jhAE2CIprf5bYDrf1ZcdUvRl9fMAwAKEIMAWncLqvwYvqXEQ+t/lrGe5XhAAaxCCAJtcOUr65ZtStCdwPDqhfpzrBAGwCBdLBGxz5Sjpe7dyxWgA1iMEATbqFCYl3xTqKgAgpPg4DAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFgpZCFo9+7dGjt2rJKTkxUREaH+/ftr+vTpqqmpCZjncDgaPV5++eUQVQ0AADqKzqF64S+//FJ1dXV65ZVXdMUVV2jz5s0aN26cTp06pblz5wbMLSwsVFZWln85JiamtcsFAAAdTMhCUFZWVkCwufzyy7Vt2za99NJLjUJQ9+7d5Xa7W7tEAADQgYUsBDWlsrJSPXv2bDSel5en+++/X8nJyRo7dqzGjx+vTp3O/0me1+uV1+v1L1dVVUmSfD6ffD5fyxfeRjX0alPPEn3Ttx3om75tEOx+HcYYE9RXaKa///3vuvbaazVv3jzdf//9/vGZM2fqlltuUUREhD766CM9/fTTmjZtmp588snzbis/P18FBQWNxhcvXqzIyMig1A8AAFpWdXW1srOzVVlZqejo6BbffouHoPMFkHMVFxdr8ODB/uWDBw8qPT1d6enpeu211y743Hnz5mnGjBmqrKw875ymjgQlJiaqrKxMsbGxzeyk/fP5fCoqKlJmZqacTmeoy2k19E3fNqBv+rZBRUWFPB5P0EJQi38clpeXpzFjxlxwTr9+/fw/Hzx4UBkZGRoyZIheffXVb91+amqqqqqqdOjQIcXHxzc5x+VyyeVyNRp3Op1WvXka0Ldd6Nsu9G0X2/oOdq8tHoLi4uIUFxfXrLkHDhxQRkaGBg0apMLCwgue59Ng06ZN6tKli7p3736JlQIAAJuF7MTogwcPaujQoerbt6/mzp2rf/zjH/51Dd8EW758ucrLyzVkyBBFRERo9erVeuKJJzR+/Pgmj/QAAAA0V8hC0MqVK7Vz507t3LlTffr0CVjXcJqS0+nUwoULNXnyZNXV1enyyy/XjBkzNHHixFCUDAAAOpCQhaDc3Fzl5uZecM43ryUEAADQUrh3GAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsFNIQ1K9fPzkcjoDH448/HjBn7969GjlypLp27aq4uDg98sgjqqmpCVHFAACgo+gc6gJmzJihcePG+Ze7devm/7m2tla33nqrevXqpXXr1qmiokI5OTkyxmjBggWhKBcAAHQQIQ9BUVFRcrvdTa5buXKltm7dqn379ikhIUGSNG/ePOXm5mrWrFmKjo5uzVIBAEAHEvIQ9Nxzz+mZZ55RYmKiRo8erd/85jcKDw+XJG3YsEEpKSn+ACRJw4cPl9frVUlJiTIyMprcptfrldfr9S9XVVVJknw+n3w+XxC7aVsaerWpZ4m+6dsO9E3fNgh2vyENQY8++qiuvfZa9ejRQ5999pmmTZumXbt26bXXXpMklZeXKz4+PuA5PXr0UHh4uMrLy8+73Tlz5qigoKDR+OrVqxUZGdmyTbQDRUVFoS4hJOjbLvRtF/q2Q3V1dVC37zDGmJbcYH5+fpMB5FzFxcUaPHhwo/H//M//1C9+8QsdOXJEsbGxGj9+vPbs2aMVK1YEzAsPD9ebb76pMWPGNLn9po4EJSYmqqysTLGxsd+hq/bJ5/OpqKhImZmZcjqdoS6n1dA3fduAvunbBhUVFfJ4PKqsrAzKKTAtfiQoLy/vvOGkQb9+/ZocT01NlSTt3LlTsbGxcrvd+vTTTwPmHDt2TD6fr9ERonO5XC65XK5G406n06o3TwP6tgt924W+7WJb38HutcVDUFxcnOLi4r7Tczdt2iRJ8ng8kqQhQ4Zo1qxZKisr84+tXLlSLpdLgwYNapmCAQCAlUJ2TtCGDRu0ceNGZWRkKCYmRsXFxfr1r3+tUaNGqW/fvpKkYcOG6corr9Q999yjf/mXf9HRo0c1depUjRs3jm+GAQCASxKyEORyubR06VIVFBTI6/UqKSlJ48aN02OPPeafExYWpvfff18TJkzQjTfeqIiICGVnZ2vu3LmhKhsAAHQQIQtB1157rTZu3Pit8/r27as///nPrVARAACwCfcOAwAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVQhaC1qxZI4fD0eSjuLjYP6+p9S+//HKoygYAAB1E51C9cFpamsrKygLGnnrqKa1atUqDBw8OGC8sLFRWVpZ/OSYmplVqBAAAHVfIQlB4eLjcbrd/2efzadmyZcrLy5PD4QiY271794C5AAAAlypkIeibli1bpiNHjig3N7fRury8PN1///1KTk7W2LFjNX78eHXqdP5P8rxer7xer3+5qqpKUn3Q8vl8LV57W9XQq009S/RN33agb/q2QbD7dRhjTFBfoZl++tOfSpI++OCDgPGZM2fqlltuUUREhD766CM9/fTTmjZtmp588snzbis/P18FBQWNxhcvXqzIyMiWLRwAAARFdXW1srOzVVlZqejo6BbffouHoPMFkHMVFxcHnPezf/9+JSUl6e2339bPf/7zCz533rx5mjFjhiorK887p6kjQYmJiSorK1NsbGwzO2n/fD6fioqKlJmZKafTGepyWg1907cN6Ju+bVBRUSGPxxO0ENTiH4fl5eVpzJgxF5zTr1+/gOXCwkLFxsZq1KhR37r91NRUVVVV6dChQ4qPj29yjsvlksvlajTudDqtevM0oG+70Ldd6NsutvUd7F5bPATFxcUpLi6u2fONMSosLNS9997brGY3bdqkLl26qHv37pdQJQAAsF3IT4z+61//ql27dmns2LGN1i1fvlzl5eUaMmSIIiIitHr1aj3xxBMaP358k0d6AAAAmivkIWjRokVKS0vT97///UbrnE6nFi5cqMmTJ6uurk6XX365ZsyYoYkTJ4agUgAA0JGEPAQtXrz4vOuysrICLpIIAADQUrh3GAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWCmoImjVrltLS0hQZGanu3bs3OWfv3r0aOXKkunbtqri4OD3yyCOqqakJmFNaWqr09HRFRETosssu04wZM2SMCWbpAACgg+sczI3X1NRo9OjRGjJkiBYtWtRofW1trW699Vb16tVL69atU0VFhXJycmSM0YIFCyRJVVVVyszMVEZGhoqLi7V9+3bl5uaqa9eumjJlSjDLBwAAHVhQQ1BBQYEk6fXXX29y/cqVK7V161bt27dPCQkJkqR58+YpNzdXs2bNUnR0tP793/9dZ86c0euvvy6Xy6WUlBRt375d8+fP1+TJk+VwOBpt1+v1yuv1+pcrKyslSUePHm3hDts2n8+n6upqVVRUyOl0hrqcVkPf9G0D+qZvGzT8ux2sT3+CGoK+zYYNG5SSkuIPQJI0fPhweb1elZSUKCMjQxs2bFB6erpcLlfAnGnTpmn37t1KTk5utN05c+b4A9i5Bg4cGJxGAABA0FRUVCgmJqbFtxvSEFReXq74+PiAsR49eig8PFzl5eX+Of369QuY0/Cc8vLyJkPQtGnTNHnyZP/y8ePHlZSUpL179wblD7GtqqqqUmJiovbt26fo6OhQl9Nq6Ju+bUDf9G2DyspK9e3bVz179gzK9i86BOXn5zd5lOVcxcXFGjx4cLO219THWcaYgPFvzmk4LNbUcyXJ5XIFHDlqEBMTY9Wbp0F0dDR9W4S+7ULfdrG1706dgvM9rosOQXl5eRozZswF53zzyM35uN1uffrppwFjx44dk8/n8x/tcbvd/qNCDQ4fPixJjY4iAQAANNdFh6C4uDjFxcW1yIsPGTJEs2bNUllZmTwej6T6k6VdLpcGDRrkn/O73/1ONTU1Cg8P989JSEhodtgCAAD4pqBeJ2jv3r36/PPPtXfvXtXW1urzzz/X559/rpMnT0qShg0bpiuvvFL33HOPNm3apI8++khTp07VuHHj/If7srOz5XK5lJubq82bN+vdd9/V7Nmzz/vNsKa4XC5Nnz69yY/IOjL6pm8b0Dd924C+g9O3wwTxqoO5ubl64403Go2vXr1aQ4cOlVQflCZMmKC//vWvioiIUHZ2tubOnRvQcGlpqSZOnKjPPvtMPXr00IMPPqinn3662SEIAADgm4IaggAAANoq7h0GAACsRAgCAABWIgQBAAArEYIAAICVOlQImjVrltLS0hQZGanu3bs3OWfv3r0aOXKkunbtqri4OD3yyCOqqakJmFNaWqr09HRFRETosssu04wZM4J287ZgWLNmjRwOR5OP4uJi/7ym1r/88sshrPzS9evXr1FPjz/+eMCc5rwH2pPdu3dr7NixSk5OVkREhPr376/p06c36qkj7u+FCxcqOTlZXbp00aBBg/TJJ5+EuqQWNWfOHF133XWKiopS7969dfvtt2vbtm0Bc3Jzcxvt19TU1BBV3DLy8/Mb9eR2u/3rjTHKz89XQkKCIiIiNHToUG3ZsiWEFbeMpv7+cjgcmjhxoqSOs68//vhjjRw5UgkJCXI4HHrvvfcC1jdn/3q9Xj388MOKi4tT165dNWrUKO3fv/+iawnpvcNaWk1NjUaPHq0hQ4Zo0aJFjdbX1tbq1ltvVa9evbRu3TpVVFQoJydHxhgtWLBAUv39WTIzM5WRkaHi4mJt375dubm56tq1q6ZMmdLaLX0naWlpKisrCxh76qmntGrVqka3MyksLFRWVpZ/uSPcW23GjBkaN26cf7lbt27+n5vzHmhvvvzyS9XV1emVV17RFVdcoc2bN2vcuHE6deqU5s6dGzC3I+3vpUuXatKkSVq4cKFuvPFGvfLKKxoxYoS2bt2qvn37hrq8FrF27VpNnDhR1113nc6ePasnnnhCw4YN09atW9W1a1f/vKysLBUWFvqXGy4s25794Ac/0KpVq/zLYWFh/p+ff/55zZ8/X6+//roGDhyomTNnKjMzU9u2bVNUVFQoym0RxcXFqq2t9S9v3rxZmZmZGj16tH+sI+zrU6dO6eqrr9Y///M/6+c//3mj9c3Zv5MmTdLy5cu1ZMkSxcbGasqUKbrttttUUlIS8F75VqYDKiwsNDExMY3GP/jgA9OpUydz4MAB/9h//Md/GJfLZSorK40xxixcuNDExMSYM2fO+OfMmTPHJCQkmLq6uqDXHgw1NTWmd+/eZsaMGQHjksy7774bmqKCJCkpyfz+978/7/rmvAc6gueff94kJycHjHW0/X399debBx98MGDse9/7nnn88cdDVFHwHT582Egya9eu9Y/l5OSYn/3sZ6ErKgimT59urr766ibX1dXVGbfbbZ599ln/2JkzZ0xMTIx5+eWXW6nC1vHoo4+a/v37+//t6Yj7+pt/LzVn/x4/ftw4nU6zZMkS/5wDBw6YTp06mQ8//PCiXr9DfRz2bTZs2KCUlBQlJCT4x4YPHy6v16uSkhL/nPT09ICLNQ4fPlwHDx7U7t27W7vkFrFs2TIdOXJEubm5jdbl5eUpLi5O1113nV5++WXV1dW1foEt7LnnnlNsbKyuueYazZo1K+Bjoea8BzqCysrKJu+63FH2d01NjUpKSjRs2LCA8WHDhmn9+vUhqir4KisrJanRvl2zZo169+6tgQMHaty4cf77K7ZnO3bsUEJCgpKTkzVmzBh99dVXkqRdu3apvLw8YN+7XC6lp6d3qH1fU1Ojt956S/fdd1/AhYE74r4+V3P2b0lJiXw+X8CchIQEpaSkXPR7oEN9HPZtysvLG910tUePHgoPD/ffpLW8vLzRPckanlNeXq7k5ORWqbUlLVq0SMOHD1diYmLA+DPPPKNbbrlFERER+uijjzRlyhQdOXJETz75ZIgqvXSPPvqorr32WvXo0UOfffaZpk2bpl27dum1116T1Lz3QHv397//XQsWLNC8efMCxjvS/j5y5Ihqa2sb7cv4+PgOsx+/yRijyZMn60c/+pFSUlL84yNGjNDo0aOVlJSkXbt26amnntLNN9+skpKSdnuLhRtuuEFvvvmmBg4cqEOHDmnmzJlKS0vTli1b/Pu3qX2/Z8+eUJQbFO+9956OHz8e8MtrR9zX39Sc/VteXq7w8HD16NGj0ZyL/f+/zYeg/Px8FRQUXHBOcXFxo3NdzqepW20YYwLGvznH/O9J0aG+Tcd3+bPYv3+/VqxYobfffrvR3HP/8bvmmmsk1Z9P09b+UbyYvn/961/7x374wx+qR48e+sUvfuE/OiQ17z3QFnyX/X3w4EFlZWVp9OjRuv/++wPmtpf9fTGa+n+1re3HlpKXl6cvvvhC69atCxi/6667/D+npKRo8ODBSkpK0vvvv68777yztctsESNGjPD/fNVVV2nIkCHq37+/3njjDf+JwB193y9atEgjRowIOGrdEff1+XyX/ftd3gNtPgTl5eVpzJgxF5zT3LvJu91uffrppwFjx44dk8/n86dOt9vdKEk2HG78ZjJtbd/lz6KwsFCxsbEaNWrUt24/NTVVVVVVOnToUMh7PdelvAca/sLcuXOnYmNjm/UeaCsutu+DBw8qIyNDQ4YM0auvvvqt22+r+7s54uLiFBYW1uT/q+2tl+Z4+OGHtWzZMn388cfq06fPBed6PB4lJSVpx44drVRd8HXt2lVXXXWVduzYodtvv11S/dEAj8fjn9OR9v2ePXu0atUqvfPOOxec1xH3dcO3AC+0f91ut2pqanTs2LGAo0GHDx9WWlraxb3gdzuVqW37thOjDx486B9bsmRJoxOju3fvbrxer3/Os88+2y5PjK6rqzPJyclmypQpzZq/YMEC06VLl4CTwtu75cuXG0lmz549xpjmvQfao/3795sBAwaYMWPGmLNnzzbrOe19f19//fXmoYceChj7/ve/36FOjK6rqzMTJ040CQkJZvv27c16zpEjR4zL5TJvvPFGkKtrPWfOnDGXXXaZKSgo8J84+9xzz/nXe73eDnVi9PTp043b7TY+n++C8zrCvtZ5Toy+0P5tODF66dKl/jkHDx78TidGd6gQtGfPHrNp0yZTUFBgunXrZjZt2mQ2bdpkTpw4YYwx5uzZsyYlJcXccsst5m9/+5tZtWqV6dOnj8nLy/Nv4/jx4yY+Pt7cfffdprS01LzzzjsmOjrazJ07N1RtfWerVq0ykszWrVsbrVu2bJl59dVXTWlpqdm5c6f5t3/7NxMdHW0eeeSREFTaMtavX2/mz59vNm3aZL766iuzdOlSk5CQYEaNGuWf05z3QHtz4MABc8UVV5ibb77Z7N+/35SVlfkfDTri/l6yZIlxOp1m0aJFZuvWrWbSpEmma9euZvfu3aEurcU89NBDJiYmxqxZsyZgv1ZXVxtjjDlx4oSZMmWKWb9+vdm1a5dZvXq1GTJkiLnssstMVVVViKv/7qZMmWLWrFljvvrqK7Nx40Zz2223maioKP++ffbZZ01MTIx55513TGlpqbn77ruNx+Np1z03qK2tNX379jW//e1vA8Y70r4+ceKE/99nSf6/txt+WW3O/n3wwQdNnz59zKpVq8zf/vY3c/PNN5urr7662b8ENuhQISgnJ8dIavRYvXq1f86ePXvMrbfeaiIiIkzPnj1NXl5eo9+Ev/jiC3PTTTcZl8tl3G63yc/Pb3dHgYwx5u677zZpaWlNrvvLX/5irrnmGtOtWzcTGRlpUlJSzAsvvPCtv3m0ZSUlJeaGG24wMTExpkuXLuaf/umfzPTp082pU6cC5jXnPdCeFBYWNvm+P/dAb0fc38YY86//+q8mKSnJhIeHm2uvvTbgq+Mdwfn2a2FhoTHGmOrqajNs2DDTq1cv43Q6Td++fU1OTo7Zu3dvaAu/RHfddZfxeDzG6XSahIQEc+edd5otW7b419fV1fmPlrhcLvPjH//YlJaWhrDilrNixQojyWzbti1gvCPt69WrVzf5vs7JyTHGNG//nj592uTl5ZmePXuaiIgIc9ttt32nPwuHMe3oUsgAAAAtxKrrBAEAADQgBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlf4/FMZ7rs/UpwYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sympy as sy\n",
    "\n",
    "d=sy.symbols('d')\n",
    "\n",
    "def Roderiguez(v1,v2):\n",
    "    \n",
    "    v1n=v1.normalized()\n",
    "    v2n=v2.normalized()\n",
    "    \n",
    "    k=v1n.cross(v2n).normalized()\n",
    "    ct=(v1n.transpose()*v2n)[0]\n",
    "    \n",
    "    st=sy.sqrt(1-ct*ct)\n",
    "    kx=k[0]\n",
    "    ky=k[1]\n",
    "    kz=k[2]\n",
    "    \n",
    "    k_as=sy.Matrix(3,3,[0,-kz,ky,kz,0,-kx,-ky,kx,0])\n",
    "    \n",
    "    I=sy.Matrix(3,3,np.eye(3).flatten())\n",
    "    \n",
    "    R=I*ct+(1-ct)*k*(k.transpose())+st*k_as\n",
    "    \n",
    "    \n",
    "    return R\n",
    "    \n",
    "v1=vc.subs('theta',0)\n",
    "\n",
    "x=np.array([[-2.01,-1.001]])*10\n",
    "\n",
    "sub_key=['p_x','p_y','p_z','p_T_x','p_T_y','p_T_z']\n",
    "sub_val=[x[0,0],x[0,1],0,0,0,0]\n",
    "sub_dict=dict(zip(sub_key,sub_val))\n",
    "\n",
    "v2=vtb.subs(sub_dict)\n",
    "R=Roderiguez(v1,v2)\n",
    "\n",
    "s_x,s_y,s_z=sy.symbols(['sx','sy','sz'])\n",
    "S=sy.Matrix(3,3,np.diag([s_x,s_y,s_z]).flatten())\n",
    "\n",
    "sub_key=['sx','sy','sz']\n",
    "sub_val=[1,0.1,2]\n",
    "sub_dict=dict(zip(sub_key,sub_val))\n",
    "\n",
    "P=(R*S*R.transpose()).subs(sub_dict)\n",
    "P_np=np.array([[P[0,0],P[0,1]],[P[1,0],P[1,1]]],dtype=np.double)\n",
    "plot_covariance_ellipse(np.zeros([2,1]),P_np)\n",
    "plt.scatter(x[0,0],x[0,1])\n",
    "#print(\"$\"+sy.latex(P[0,0])+\"$\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 2)\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "index 1 is out of bounds for axis 0 with size 1",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_21756/1712742207.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     35\u001b[0m \u001b[0mxEst\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\u001b[0m \u001b[0mPEst\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdiag\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m40\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 37\u001b[0;31m \u001b[0mplot_covariance_ellipse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mPEst\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/tmp/ipykernel_21756/1712742207.py\u001b[0m in \u001b[0;36mplot_covariance_ellipse\u001b[0;34m(xEst, PEst)\u001b[0m\n\u001b[1;32m     25\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m     \u001b[0mpx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mxEst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 27\u001b[0;31m     \u001b[0mpy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mxEst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     28\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxEst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mxEst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     29\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"--r\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1 is out of bounds for axis 0 with size 1"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from scipy.spatial.transform import Rotation as Rot\n",
    "import matplotlib.pyplot as plt\n",
    "def plot_covariance_ellipse(xEst, PEst):  # pragma: no cover\n",
    "    Pxy = PEst[0:2, 0:2]\n",
    "    eigval, eigvec = np.linalg.eig(Pxy)\n",
    "\n",
    "    if eigval[0] >= eigval[1]:\n",
    "        bigind = 0\n",
    "        smallind = 1\n",
    "    else:\n",
    "        bigind = 1\n",
    "        smallind = 0\n",
    "\n",
    "    t = np.arange(0, 2 * math.pi + 0.1, 0.1)\n",
    "    a = math.sqrt(eigval[bigind])\n",
    "    b = math.sqrt(eigval[smallind])\n",
    "    x = [a * math.cos(it) for it in t]\n",
    "    y = [b * math.sin(it) for it in t]\n",
    "    angle = math.atan2(eigvec[1, bigind], eigvec[0, bigind])\n",
    "    rot = Rot.from_euler('z', angle).as_matrix()[0:2, 0:2]\n",
    "    fx = rot @ (np.array([x, y]))*10\n",
    "    \n",
    "    print(xEst.shape)\n",
    "    \n",
    "    px = np.array(fx[0, :]+xEst[0,0]).flatten()\n",
    "    py = np.array(fx[1, :]+xEst[1,0]).flatten()\n",
    "    plt.scatter(xEst[0,0],xEst[1,0])\n",
    "    plt.plot(px, py, \"--r\")\n",
    "    plt.grid()\n",
    "    #plt.xlim([-100,100])\n",
    "    #plt.ylim([-100,100])\n",
    "    \n",
    "    \n",
    "xEst=np.ones([2,1])\n",
    "PEst=np.diag([5,40])\n",
    "plot_covariance_ellipse(xEst,PEst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.patches.FancyArrow at 0x7febcb24bf40>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3kklEQVR4nO3deXRUdYL+/ycQCNlZA0knaJAdCa0JQhwI2EEgODQobU8rKirNGb+yqIh6QG3bM2icXyOuLK3gRtpGlGaxWVrQJsGRKGERhIhxCCRCFrbsoZKqur8/aO4YJZDKdmt5v86pPz6Ve6uepCD15PO5da+fYRiGAAAALNLG6gAAAMC3UUYAAIClKCMAAMBSlBEAAGApyggAALAUZQQAAFiKMgIAACxFGQEAAJbytzpAQzidTp08eVKhoaHy8/OzOg4AAGgAwzBUXl6uqKgotWlT//yHR5SRkydPKiYmxuoYAACgEfLz8xUdHV3v1z2ijISGhkq68M2EhYVZnAYAADREWVmZYmJizPfx+nhEGbm4NBMWFkYZAQDAw1zpEAsOYAUAAJaijAAAAEtRRgAAgKUoIwAAwFKUEQAAYCnKCAAAsBRlBAAAWIoyAgAALEUZAQAAlqKMAAAAS1FGAACApSgjAADAUpQRAAB8WFVVlWw2m6UZKCMAAPiYyspKffTRR/rtb3+rrt26afZDD1uax9/SZwcAAK2isrJSmzZt0po1a7Rp8xadr65Sm/Yd5Kw5r0ED+luajTICAICXqqio0N///net+fBDbd68Rbbz1QqM7K0ON9yujlcN0dkNz+vGf7tRs2fPtjQnZQQAAC9SXl6uv//97/pgzRpt2bJVNbbzCozqq8Dh/6Eu/UaoXccekqTTG/8/BRi1WvXeu2rTxtqjNigjAAB4uLKyMn388cf6YM0abd36D9XW2BT4i34KSrxTXfvdaBaQiyoP71Bldob++te/Kjo62qLU/4cyAgCAByotLdXGjRu15sMP9Y9/fPKvAtJfIf82VUH9Rsg/POKS+9nLTqlk+3L9x+9+p9/97netnPrSKCMAAHiIkpISbdy4UR+sWaNtn2xTbW2NAqMHKGTE3Qrqd6P8wy5dQC4yDKfObX5ZXTuFa9nSpa2U+sooIwAAuLFz585pw4YN+mDNGm3ftl12e62CYgYqJGmagvreKP+wbg1+rPKsjao6/rU2bt+uTp06tWBq11BGAABwM2fPntX69eu15sMP9en2T2V32BUUPVCho+67MAMS2tXlx6w5dUylGe/poYceUnJycgukbjzKCAAAbuDMmTNav369PlizRp999pkcDoeCYgYpdPT9F2ZAQrs0+rENe63ObV6sPr2v0QsvvNCMqZsHZQQAAAsZhqFbb5uijz/eKMNpKLDntQq/6fcK7Huj/EM6N8tzlHz+F9lP52v1lq/UoUOHZnnM5kQZAQDAQuXl5frHP/6htp2j1f0/FqptcPMey3E+/xuVfbVWL6Sm6pe//GWzPnZz4do0AABYKCwsTMuXLVXtqeOqPrq3WR/baatSyeaXlJh4o+bNm9esj92cKCMAAFhs2rRpuvfee1WybalqTh1vtsc9u/3Pamuv1F/SVqlt27bN9rjNjTICAIAbWLJkifr06a1zH/+3nDXVTX68yiP/o8pvPtWS115TbGxsMyRsOU0qI6mpqfLz89PDDz982e3S09MVHx+vDh06qFevXlq+fHlTnhYAAK8TFBSkdWs/kl/lGZ39ZKkMw2j0Y9krzqp021JNmjxZ06ZNa8aULaPRZWT37t164403FBcXd9ntcnNzNWHCBI0cOVL79u3TggULNGfOHK1du7axTw0AgFfq37+/Vrz5hioP/VMVBz5p1GMYhqFzW19VeFCAVrz5pvz8/Jo5ZfNrVBmpqKjQ1KlT9eabb17xDG7Lly9Xz5499fLLL2vAgAH6/e9/r/vvv1+LFi1qVGAAALzZnXfeqRkzZqj00zdUU5zr8v4V+7eo6n+z9O47b6trV9dPjmaFRpWRmTNn6pZbbtGYMWOuuO2uXbs0duzYOveNGzdOWVlZqq2tveQ+NptNZWVldW4AAPiKV155RQP699O5jS/Iaatq8H61Z35Q6Y6V+s///E9NmDChBRM2L5fLyOrVq7V3716lpqY2aPvCwkJ17969zn3du3eX3W7X6dOnL7lPamqqwsPDzVtMTIyrMQEA8Fjt27fXiR/ydf7MCRW8M7tBx48YDrvObV6sq2Ji9OKLL7ZCyubjUhnJz8/XQw89pLS0NJfO4PbT9aqLP9T61rHmz5+v0tJS85afn+9KTAAAPFZubq78/f117tw5SZK9pEgV+7dccb/SXR/IVvS/ev8vaQoODm7pmM3KpTKyZ88eFRcXKz4+Xv7+/vL391d6erpeffVV+fv7y+Fw/GyfHj16qLCwsM59xcXF8vf3V5culz7PfkBAgMLCwurcAADwdq+99pp69eolSbrrrrvkdDr14IMPquSzN2Ur/L7e/Wwnj6hs1xo9/dRTGjZsWGvFbTYunQ4+OTlZBw8erHPffffdp/79++uJJ5645AlVEhMT9fHHH9e575NPPlFCQoLatWvXiMgAAHgXh8OhXr16KS8vT5L02Wef6aabbpIkLV68WF/sylT2x/+tiHteVpuAurMezprzOrd5sa6//no9+eSTrZ69Obg0MxIaGqprr722zi04OFhdunTRtddeK+nCEss999xj7vPAAw/o+PHjmjt3rrKzs/XWW29p5cqVbn1aWgAAWktOTo78/f3NIlJeXm4WEenCasHf1n6kdvZKnd3yys+OHzn3zxXyqzqn9/+S5rF/5Df7GVgLCgrMH6gkxcbGavPmzdqxY4d++ctf6r/+67/06quvasqUKc391AAAeJRFixapb9++kqQZM2bIMAyFhIT8bLvY2Fi99847qjzyhcr3/N9qQ9X/7lbF/q16afGL5uN4Ij+jKad4ayVlZWUKDw9XaWkpx48AADxebW2toqKizE+V7ty5UyNGjLjifg8//LBee32JIu78b/l37KHid2bppn8bpi2bN7vlyc0a+v5NGQEAoBUdPnxYgwYNMseVlZUKCgpq0L41NTW68d9G6Jvv89S2S4zanz2qw4e+UWRkZEvFbZKGvn9zoTwAAFrJwoULzSIyZ84cGYbR4CIiXTj/yNqPPlSAYVPV/2Zp5Yo33baIuMKlT9MAAADX1dTUqFOnTqqqunA21czMzEZ/BPeqq67SPz/7VCdPntS///u/N2dMy1BGAABoQQcOHNCQIUPMcXV1tUsnDr2U66+/Xtdff31To7kNlmkAAGghTz31lFlEHn/8cRmG0eQi4o2YGQEAoJmdP39egYGB5njv3r267rrrLEzk3pgZAQCgGWVlZZlFpH379jp//jxF5AooIwAANJNHH31UQ4cOlSQ988wzstlsCggIsDiV+2OZBgCAJqqurq7zEd0DBw5o8ODBFibyLMyMAADQBLt27TKLSHh4uGpqaigiLqKMAADQSDNnztSNN94oSUpNTVVJSYnHXqzOSizTAADgosrKyjoXtMvOzlb//v0tTOTZmBkBAMAF6enpZhGJiopSbW0tRaSJKCMAADTQ/fffr9GjR0uSXnrpJZ04cUL+/iwyNBU/QQAAruDi1WcvysnJUe/evS1M5F2YGQEA4DK2b99uFpFrrrlGdrudItLMKCMAAFyCYRi64447dPPNN0uSli5dqu+//15t27a1OJn3YZkGAICfKCkpUadOncxxbm6urr76ausCeTlmRgAA+JHNmzebRWTw4MFyOBwUkRZGGQEAQBeWZSZNmqRbbrlFkrRy5UodOHBAbdrwVtnSWKYBAPi8M2fOqGvXruY4Pz9f0dHRFibyLdQ9AIBPW7dunVlEhg0bJofDQRFpZZQRAIBPMgxDY8eO1W233SZJSktLU2ZmJssyFmCZBgDgc4qLi9W9e3dzXFBQoB49eliYyLdR/wAAPuWDDz4wi8jo0aPldDopIhZjZgQA4BMMw1BSUpI+//xzSdKHH36o3/zmNxangkQZAQD4gJMnT+oXv/iFOS4uLla3bt0sTIQfY5kGAODV3nvvPbOITJgwQU6nkyLiZpgZAQB4JafTqYSEBO3bt0+StGHDBv3617+2OBUuhTICAPA6+fn56tmzpzk+c+aMOnfubGEiXA7LNAAAr/Lmm2+aRWTKlClyOp0UETfHzAgAwCs4HA5de+21+vbbbyVJW7Zs0fjx4y1OhYagjAAAPF5ubq569epljktKShQeHm5hIriCZRoAgEd77bXXzCJy1113yel0UkQ8DDMjAACP5HA41KtXL+Xl5UmSPv30U/3qV7+yOBUagzICAPA4OTk56tu3rzkuLy9XSEiIhYnQFCzTAAA8yqJFi8wiMmPGDBmGQRHxcMyMAAA8Qm1traKionT69GlJ0s6dOzVixAiLU6E5UEYAAG4vOztbAwcONMeVlZUKCgqyMBGaE8s0AAC39txzz5lFZM6cOTIMgyLiZZgZAQC4pZqaGnXq1ElVVVWSpMzMTA0bNsziVGgJlBEAgNs5cOCAhgwZYo6rq6vVoUMHCxOhJbm0TLNs2TLFxcUpLCxMYWFhSkxM1JYtW+rdfseOHfLz8/vZ7eKpegEA+KmnnnrKLCKPP/64DMOgiHg5l2ZGoqOj9cILL6h3796SpHfffVeTJk3Svn37NGjQoHr3O3LkiMLCwsxxt27dGhkXAOCtbDZbndKxZ88eXX/99RYmQmtxqYxMnDixzvi5557TsmXLlJmZedkyEhERoY4dOzYqIADA++3Zs0cJCQmSJH9/f1VUVCggIMDiVGgtjf40jcPh0OrVq1VZWanExMTLbnvdddcpMjJSycnJ+uc//3nFx7bZbCorK6tzAwB4p3nz5plF5A9/+INqa2spIj7G5QNYDx48qMTERJ0/f14hISFat25dnc9+/1hkZKTeeOMNxcfHy2azadWqVUpOTtaOHTuUlJRU73Okpqbq2WefdTUaAMCDVFdX1/mI7oEDBzR48GALE8EqfoZhGK7sUFNTo7y8PJWUlGjt2rVasWKF0tPT6y0kPzVx4kT5+flp48aN9W5js9lks9nMcVlZmWJiYlRaWlrn2BMAgGfKzMw0Z9XDw8N16tQptWvXzuJUaG5lZWUKDw+/4vu3y8s07du3V+/evZWQkKDU1FQNGTJEr7zySoP3Hz58uHJyci67TUBAgPmJnYs3AIB3mDVrlllEUlNTVVJSQhHxcU0+z4hhGHVmMa5k3759ioyMbOrTAgA8TGVlZZ0L2mVnZ6t///4WJoK7cKmMLFiwQCkpKYqJiVF5eblWr16tHTt2aOvWrZKk+fPn68SJE3rvvfckSS+//LKuvvpqDRo0SDU1NUpLS9PatWu1du3a5v9OAABuKyMjQ6NGjZJ04XjCvLw8+ftz3k1c4NK/hKKiIt19990qKChQeHi44uLitHXrVt18882SpIKCAuXl5Znb19TUaN68eTpx4oQCAwM1aNAgbdq0SRMmTGje7wIA4Lbuv/9+vf3225KkxYsX65FHHrE4EdyNywewWqGhB8AAANxHeXl5nd/ZOTk55kkz4Rta7ABWAACuZPv27eabT2xsrOx2O0UE9aKMAACajWEYuvPOO83l+6VLl+ro0aNq27atxcngzjh6CADQLEpKStSpUydznJubq6uvvtq6QPAYzIwAAJps8+bNZhEZPHiwHA4HRQQNRhkBADSaYRiaPHmybrnlFknSypUrdeDAAbVpw9sLGo5lGgBAo5w5c0Zdu3Y1x/n5+YqOjrYwETwV1RUA4LL169ebReSGG26Qw+GgiKDRKCMAgAYzDEPjxo3TrbfeKklKS0vTl19+ybIMmoRlGgBAgxQXF6t79+7muKCgQD169LAwEbwFVRYAcEVr1qwxi8ioUaPkdDopImg2zIwAAOplGIaSkpL0+eefS5I+/PBD/eY3v7E4FbwNZQQAcEkFBQWKiooyx8XFxerWrZuFieCtWKYBAPzMqlWrzCKSkpIip9NJEUGLYWYEAGByOp0aOnSo9u7dK0nasGGDfv3rX1ucCt6OMgIAkHThpGU9e/Y0x2fOnFHnzp0tTARfwTINAEArVqwwi8htt90mp9NJEUGrYWYEAHyYw+HQ4MGDlZ2dLUnasmWLxo8fb3Eq+BrKCAD4qNzcXPXq1cscl5SUKDw83MJE8FUs0wCAD1qyZIlZRO666y45nU6KCCzDzAgA+BCHw6FrrrlGx48flyR9+umn+tWvfmVxKvg6yggA+IicnBz17dvXHJeVlSk0NNTCRMAFLNMAgA948cUXzSLy+9//XoZhUETgNpgZAQAvZrfbFRUVpVOnTkmSdu7cqREjRlicCqiLMgIAXio7O1sDBw40xxUVFQoODrYwEXBpLNMAgBd67rnnzCIyZ84cGYZBEYHbYmYEALxITU2NOnXqpKqqKklSZmamhg0bZnEq4PIoIwDgJQ4cOKAhQ4aY4+rqanXo0MHCREDDsEwDAF7g6aefNovI448/LsMwKCLwGMyMAIAHs9lsCgwMlGEYkqQ9e/bo+uuvtzgV4BpmRgDAQ+3Zs0cdOnSQYRjy9/fX+fPnKSLwSJQRAPBAjz32mBISEiRdWKKpra1VQECAxamAxmGZBgA8SHV1tYKCgszxgQMHNHjwYAsTAU3HzAgAeIjMzEyziISFhammpoYiAq9AGQEADzBr1iwlJiZKkp5//nmVlpaqXbt2FqcCmgfLNADgxiorKxUSEmKOs7Oz1b9/fwsTAc2PmREAcFMZGRlmEenRo4dqa2spIvBKlBEAcEPTp0/XqFGjJEmLFy9WQUGB/P2ZzIZ34l82ALiR8vJyhYWFmeOcnBz17t3bwkRAy2NmBADcxKeffmoWkdjYWNntdooIfAJlBAAsZhiG7rzzTo0ZM0aStHTpUh09elRt27a1OBnQOlimAQALlZSUqFOnTuY4NzdXV199tXWBAAswMwIAFtmyZYtZRAYNGiSHw0ERgU9yqYwsW7ZMcXFxCgsLU1hYmBITE7Vly5bL7pOenq74+Hh16NBBvXr10vLly5sUGAA8nWEYmjx5siZMmCBJWrFihb755hu1acPfh/BNLi3TREdH64UXXjAPqHr33Xc1adIk7du3T4MGDfrZ9rm5uZowYYJmzJihtLQ0/c///I8efPBBdevWTVOmTGme7wAAPMiZM2fUtWtXc5yfn6/o6GgLEwHW8zMMw2jKA3Tu3Fl/+tOfNH369J997YknntDGjRuVnZ1t3vfAAw/o66+/1q5duxr8HGVlZQoPD1dpaWmdj7wBgCdZv369br31VknS0KFDlZmZyWwIvFpD378b/b/A4XBo9erVqqysNK+X8FO7du3S2LFj69w3btw4ZWVlqba2tt7HttlsKisrq3MDAE9lGIbGjRtnFpFVq1bpq6++oogA/+Lyp2kOHjyoxMREnT9/XiEhIVq3bp0GDhx4yW0LCwvVvXv3Ovd1795ddrtdp0+fVmRk5CX3S01N1bPPPutqNABwO8XFxXV+D548ebLe332Ar3K5lvfr10/79+9XZmam/t//+3+aNm2aDh8+XO/2fn5+dcYXV4V+ev+PzZ8/X6WlpeYtPz/f1ZgAYLk1a9aYRSQpKUlOp5MiAlyCyzMj7du3Nw9gTUhI0O7du/XKK6/oz3/+88+27dGjhwoLC+vcV1xcLH9/f3Xp0qXe5wgICFBAQICr0QDALRiGoVGjRmnnzp2SLpSS22+/3eJUgPtq8knPDMOQzWa75NcSExP18ccf17nvk08+UUJCgtq1a9fUpwYAt1NQUKCoqChzXFxcrG7dulmYCHB/Li3TLFiwQDt37tSxY8d08OBBPfnkk9qxY4emTp0q6cLyyj333GNu/8ADD+j48eOaO3eusrOz9dZbb2nlypWaN29e834XAOAGVq1aZRaR8ePHy+l0UkSABnBpZqSoqEh33323CgoKFB4erri4OG3dulU333yzpAt/EeTl5Znbx8bGavPmzXrkkUe0ZMkSRUVF6dVXX+UcIwC8itPp1A033KA9e/ZIuvAR3kmTJlmcCvAcTT7PSGvgPCMA3FV+fr569uxpjs+cOaPOnTtbmAhwHy1+nhEA8HUrVqwwi8htt90mp9NJEQEagav2AoCLHA6HBg8ebJ5desuWLRo/frzFqQDPRRkBABfk5uaqV69e5rikpETh4eEWJgI8H8s0ANBAS5YsMYvI1KlT5XQ6KSJAM2BmBACuwOFwqHfv3jp27Jgkafv27UpOTrY2FOBFKCMAcBk5OTnq27evOS4rK1NoaKiFiQDvwzINANTjxRdfNIvI9OnTZRgGRQRoAcyMAMBP2O12RUVF6dSpU5KkjIwMjRw50uJUgPeijADAj2RnZ2vgwIHmuKKiQsHBwRYmArwfyzQA8C/PPfecWURmzZolwzAoIkArYGYEgM+rqalR586dVVlZKUnKzMzUsGHDLE4F+A7KCACfduDAAQ0ZMsQcV1VVKTAw0MJEgO9hmQaAz3r66afNIvLYY4/JMAyKCGABZkYA+BybzabAwEBdvGj5nj17dP3111ucCvBdzIwA8Cl79uxRhw4dZBiG2rZtq/Pnz1NEAItRRgD4jMcee0wJCQmSLizR2O12BQQEWJwKAMs0ALxedXW1goKCzPGBAwc0ePBgCxMB+DFmRgB4tczMTLOIhIWFqaamhiICuBnKCACvNXv2bCUmJkq6cEKz0tJStWvXzuJUAH6KZRoAXqeyslIhISHm+PDhwxowYICFiQBcDjMjALzKzp07zSLSo0cP1dbWUkQAN0cZAeA1pk+frqSkJEnSiy++qIKCAvn7MwEMuDv+lwLweOXl5QoLCzPH3333nfr06WNhIgCuYGYEgEf79NNPzSISGxsru91OEQE8DGUEgEcyDEN33nmnxowZI0lasmSJjh49qrZt21qcDICrWKYB4HFKSkrUqVMnc5ybm6urr77aukAAmoSZEQAeZcuWLWYRGTRokBwOB0UE8HCUEQAewTAM3XrrrZowYYIkacWKFfrmm2/Upg2/xgBPxzINALd35swZde3a1Rzn5+crOjrawkQAmhN/UgBwaxs2bDCLSEJCghwOB0UE8DKUEQBuyTAMjR8/XpMnT5YkrVq1Srt372ZZBvBCLNMAcDvFxcXq3r27OT558qQiIyMtTASgJfEnBgC3smbNGrOIJCUlyel0UkQAL8fMCAC3YBiGRo0apZ07d0q6UEpuv/12i1MBaA2UEQCWKygoUFRUlDkuKipSRESEhYkAtCaWaQBYKi0tzSwi48ePl9PppIgAPoaZEQCWcDqduuGGG7Rnzx5J0vr16zVp0iSLUwGwAmUEQKvLz89Xz549zfGZM2fUuXNnCxMBsBLLNABa1YoVK8wicuutt8rpdFJEAB/HzAiAVuFwOBQXF6fDhw9LkjZv3qyUlBSLUwFwB5QRAC0uNzdXvXr1MsclJSUKDw+3MBEAd8IyDYAWtWTJErOITJ06VU6nkyICoA6XykhqaqqGDh2q0NBQRUREaPLkyTpy5Mhl99mxY4f8/Px+dvv222+bFByAe3M4HIqNjdWsWbMkSdu3b1daWpr8/PwsTgbA3bi0TJOenq6ZM2dq6NChstvtevLJJzV27FgdPnxYwcHBl933yJEjCgsLM8fdunVrXGIAbi8nJ0d9+/Y1x2VlZQoNDbUwEQB35lIZ2bp1a53x22+/rYiICO3Zs0dJSUmX3TciIkIdO3Z0OSAAz7J48WI9+uijkqTp06drxYoVFicC4O6adABraWmpJDXoY3nXXXedzp8/r4EDB+qpp57STTfdVO+2NptNNpvNHJeVlTUlJoBWYLfb9Ytf/ELFxcWSpIyMDI0cOdLiVAA8QaMPYDUMQ3PnztWIESN07bXX1rtdZGSk3njjDa1du1Z/+9vf1K9fPyUnJysjI6PefVJTUxUeHm7eYmJiGhsTQCvIzs5Wu3btzCJSUVFBEQHQYH6GYRiN2XHmzJnatGmTPv/8c0VHR7u078SJE+Xn56eNGzde8uuXmhmJiYlRaWlpneNOAFjv+eef15NPPilJmjVrll577TWLEwFwF2VlZQoPD7/i+3ejlmlmz56tjRs3KiMjw+UiIknDhw9XWlpavV8PCAhQQEBAY6IBaCW1tbXq3LmzKioqJEm7du3S8OHDLU4FwBO5VEYMw9Ds2bO1bt067dixQ7GxsY160n379ikyMrJR+wKw3sGDBxUXF2eOq6qqFBgYaGEiAJ7MpWNGZs6cqbS0NL3//vsKDQ1VYWGhCgsLVV1dbW4zf/583XPPPeb45Zdf1vr165WTk6NDhw5p/vz5Wrt2rXnuAQCe5Q9/+INZRB577DEZhkERAdAkLs2MLFu2TJI0evToOve//fbbuvfeeyVJBQUFysvLM79WU1OjefPm6cSJEwoMDNSgQYO0adMmTZgwoWnJAbQqm82moKAgOZ1OSVJWVpbi4+MtTgXAGzT6ANbW1NADYAC0jD179ighIUGS1LZtW1VWVnJcF4Arauj7N9emAXBZjz32mFlEnn76adntdooIgGbFVXsBXFJ1dbWCgoLM8ddff13noFUAaC7MjAD4mS+//NIsIqGhobLZbBQRAC2GMgKgjtmzZ5vnC3nuuedUVlam9u3bW5wKgDdjmQaAJKmyslIhISHm+PDhwxowYICFiQD4CmZGAGjnzp1mEYmIiFBtbS1FBECroYwAPm769OlKSkqSJL344osqKiqSvz+TpgBaD79xAB9VXl5e53P/3333nfr06WNhIgC+ipkRwAd9+umnZhGJjY2V3W6niACwDGUE8CGGYWjq1KkaM2aMJOn111/X0aNH1bZtW4uTAfBlLNMAPqKkpESdOnUyx7m5ubr66qutCwQA/8LMCOADtmzZYhaRgQMHym63U0QAuA3KCODFDMPQbbfdZl4le8WKFTp06BDLMgDcCss0gJc6e/asunTpYo7z8vIUExNjYSIAuDRmRgAvtGHDBrOIJCQkyOFwUEQAuC3KCOBFDMPQ+PHjNXnyZEnSqlWrtHv3brVpw391AO6LZRrAS5w6dUoRERHm+OTJk4qMjLQwEQA0DH8uAV7gww8/NItIUlKSnE4nRQSAx2BmBPBghmFo9OjRysjIkCR98MEH+u1vf2txKgBwDWUE8FAFBQWKiooyx0VFRXWWaQDAU7BMA3igtLQ0s4iMGzdOTqeTIgLAYzEzAngQp9OpYcOGKSsrS5K0fv16TZo0yeJUANA0lBHAQ/zwww91zhVy+vTpOic1AwBPxTIN4AFWrlxpFpFbb71VTqeTIgLAazAzArgxh8OhuLg4HT58WJK0efNmpaSkWJwKAJoXZQRwU8eOHVNsbKw5PnfunDp27GhdIABoISzTAG5oyZIlZhG544475HQ6KSIAvBYzI4AbcTgc6tOnj3JzcyVJ27dvV3JyssWpAKBlUUYAN/H999+rT58+5risrEyhoaEWJgKA1sEyDeAGFi9ebBaR6dOnyzAMiggAn8HMCGAhu92uX/ziFyouLpYkZWRkaOTIkRanAoDWRRkBLJKdna2BAwea44qKCgUHB1uYCACswTINYIHnn3/eLCIzZ86UYRgUEQA+i5kRoBXV1taqc+fOqqiokCTt2rVLw4cPtzgVAFiLMgK0koMHDyouLs4cV1VVKTAw0MJEAOAeWKYBWsEf/vAHs4jMmzdPhmFQRADgX5gZAVqQzWZTUFCQnE6nJCkrK0vx8fEWpwIA90IZAVrI3r17zeLRpk0bVVVVKSAgwOJUAOB+WKYBWsDjjz9uFpGnn35aDoeDIgIA9WBmBGhG1dXVCgoKMsdff/11nYNWAQA/x8wI0Ey+/PJLs4iEhITIZrNRRACgASgjQDOYPXu2eb6QhQsXqry8XO3bt7c4FQB4BpfKSGpqqoYOHarQ0FBFRERo8uTJOnLkyBX3S09PV3x8vDp06KBevXpp+fLljQ4MuJPKykr5+fnp9ddflyQdPnxYTz75pMWpAMCzuFRG0tPTNXPmTGVmZmrbtm2y2+0aO3asKisr690nNzdXEyZM0MiRI7Vv3z4tWLBAc+bM0dq1a5scHrDSzp07FRISIkmKiIhQbW2tBgwYYHEqAPA8foZhGI3d+dSpU4qIiFB6erqSkpIuuc0TTzyhjRs3Kjs727zvgQce0Ndff61du3Y16HnKysoUHh6u0tJShYWFNTYu0GxmzJihFStWSJJefPFFzZ071+JEAOB+Gvr+3aRP05SWlkqSOnfuXO82u3bt0tixY+vcN27cOK1cuVK1tbVq167dz/ax2Wyy2WzmuKysrCkxgWZTXl5e5z/Ud999pz59+liYCAA8X6MPYDUMQ3PnztWIESN07bXX1rtdYWGhunfvXue+7t27y2636/Tp05fcJzU1VeHh4eYtJiamsTGBZvPZZ5+ZReSqq66S3W6niABAM2h0GZk1a5YOHDigv/71r1fc1s/Pr8744srQT++/aP78+SotLTVv+fn5jY0JNJlhGJo6daqSk5MlSa+//rqOHTumtm3bWpwMALxDo5ZpZs+erY0bNyojI0PR0dGX3bZHjx4qLCysc19xcbH8/f3VpUuXS+4TEBDA2SrhFkpLS9WxY0dzfPToUcXGxloXCAC8kEszI4ZhaNasWfrb3/6mzz77rEG/lBMTE7Vt27Y6933yySdKSEi45PEigLvYunWrWUQGDBggu91OEQGAFuBSGZk5c6bS0tL0/vvvKzQ0VIWFhSosLFR1dbW5zfz583XPPfeY4wceeEDHjx/X3LlzlZ2drbfeeksrV67UvHnzmu+7AJqRYRi67bbblJKSIkl68803dfjwYZZlAKCFuLRMs2zZMknS6NGj69z/9ttv695775UkFRQUKC8vz/xabGysNm/erEceeURLlixRVFSUXn31VU2ZMqVpyYEWcPbs2TrLh3l5eRxADQAtrEnnGWktnGcErWHDhg2aPHmyJCk+Pl5fffWV2rThigkA0FgNff/mNy18nmEYGj9+vFlE3nvvPWVlZVFEAKCVNOmkZ4Cnu3gW4YtOnjypyMhICxMBgO/hTz/4rA8//NAsIiNHjpTT6aSIAIAFKCPwOYZhaNSoUfrtb38rSVq9erUyMjLqPQkfAKBlsUwDn1JYWFhn9qOoqKjOMg0AoPUxMwKfkZaWZhaRcePGyel0UkQAwA0wMwKv53Q6NXz4cO3evVuStH79ek2aNMniVACAiygj8Go//PBDnZOWnT59ut5rIgEArMEyDbzWypUrzSIyadIkOZ1OiggAuCFmRuB1HA6HhgwZokOHDkmSNm/ebF5nBgDgfigj8CrHjh2rc2Xdc+fOmVfeBQC4J5Zp4DWWLl1qFpE77rhDTqeTIgIAHoCZEXg8h8OhPn36KDc3V5K0fft2JScnW5wKANBQlBF4tO+//159+vQxx2VlZQoNDbUwEQDAVSzTwGMtXrzYLCL333+/DMOgiACAB2JmBB7HbrcrOjpaRUVFkqT09HQlJSVZnAoA0FiUEXiUb7/9VgMGDDDHFRUVCg4OtjARAKCpWKaBx0hNTTWLyMyZM2UYBkUEALwAMyNwe7W1terSpYvKy8slSV988YUSExMtTgUAaC6UEbi1gwcPKi4uzhxXVVUpMDDQwkQAgObGMg3c1jPPPGMWkUcffVSGYVBEAMALMTMCt2Oz2RQcHCyHwyFJysrKUnx8vMWpAAAthTICt7J3716zeLRp00ZVVVUKCAiwOBUAoCWxTAO38fjjj5tF5KmnnpLD4aCIAIAPYGYElquurlZQUJA5/vrrr+sctAoA8G7MjMBSX375pVlEQkJCZLPZKCIA4GMoI7DMnDlzNHz4cEnSwoULVV5ervbt21ucCgDQ2limQaurrKxUSEiIOT58+HCdU7wDAHwLMyNoVTt37jSLSEREhGpraykiAODjKCNoNTNmzDCvrrto0SIVFRXJ35/JOQDwdbwToMWVl5crLCzMHH/33Xfq06ePhYkAAO6EmRG0qM8++8wsIldddZXsdjtFBABQB2UELcIwDN11111KTk6WJL322ms6duyY2rZta3EyAIC7YZkGza60tFQdO3Y0x0ePHlVsbKx1gQAAbo2ZETSrrVu3mkVkwIABstvtFBEAwGVRRtAsDMPQbbfdppSUFEnSm2++qcOHD7MsAwC4IpZp0GRnz55Vly5dzHFeXp5iYmIsTAQA8CTMjKBJNm7caBaR+Ph4ORwOiggAwCWUETSKYRhKSUnRpEmTJEnvvvuusrKy1KYN/6QAAK5hmQYuO3XqlCIiIszxiRMnFBUVZWEiAIAn489YuOSjjz4yi8iIESPkdDopIgCAJnG5jGRkZGjixImKioqSn5+f1q9ff9ntd+zYIT8/v5/dvv3228ZmhgUMw9CoUaN0++23S5JWr16tnTt3ys/Pz+JkAABP5/IyTWVlpYYMGaL77rtPU6ZMafB+R44cqXN9km7durn61LBIYWGhIiMjzXFRUVGdZRoAAJrC5TKSkpJinkvCFREREXXOygnPkJaWprvvvluSNHbsWG3dupXZEABAs2q1Y0auu+46RUZGKjk5Wf/85z9b62nRSE6nUzfccINZRNatW6d//OMfFBEAQLNr8U/TREZG6o033lB8fLxsNptWrVql5ORk7dixQ0lJSZfcx2azyWazmeOysrKWjokf+eGHH+qcK+T06dN1TmoGAEBzavEy0q9fP/Xr188cJyYmKj8/X4sWLaq3jKSmpurZZ59t6Wi4hLfeekvTp0+XJE2aNEnr1q1jNgQA0KIs+Wjv8OHDlZOTU+/X58+fr9LSUvOWn5/fiul8k9Pp1ODBg80ismnTJq1fv54iAgBocZac9Gzfvn11Pp3xUwEBAQoICGjFRL7t2LFjda6se+7cOQ42BgC0GpdnRioqKrR//37t379fkpSbm6v9+/crLy9P0oVZjXvuucfc/uWXX9b69euVk5OjQ4cOaf78+Vq7dq1mzZrVPN8BmmTp0qVmEbnjjjvkdDopIgCAVuXyzEhWVpZuuukmczx37lxJ0rRp0/TOO++ooKDALCaSVFNTo3nz5unEiRMKDAzUoEGDtGnTJk2YMKEZ4qOxHA6H+vTpo9zcXEnStm3bNGbMGItTAQB8kZ9hGIbVIa6krKxM4eHhKi0trXPiNDTO999/rz59+phjfq4AgJbQ0Pdvrk3jY1566SWziNx3330yDIMiAgCwFFft9RF2u13R0dEqKiqSJKWnp9f70WoAAFoTZcQHfPvttxowYIA5rqioUHBwsIWJAAD4PyzTeLnU1FSziDz44IMyDIMiAgBwK8yMeKna2lp17drVPJX+F198ocTERItTAQDwc5QRL3Tw4EHFxcWZ46qqKgUGBlqYCACA+rFM42WeeeYZs4g8+uijMgyDIgIAcGvMjHgJm82m4OBgORwOSRdOThcfH29xKgAArowy4gX27t1rFg8/Pz9VV1dzbR8AgMdgmcbDPfHEE2YRefLJJ+V0OikiAACPwsyIh6qurlZQUJA53r9/v4YMGWJhIgAAGoeZEQ/05ZdfmkUkODhYNpuNIgIA8FiUEQ8zZ84cDR8+XJK0cOFCVVRUqH379hanAgCg8Vim8RCVlZUKCQkxx4cOHdLAgQMtTAQAQPNgZsQDfP7552YR6datm2praykiAACvQRlxczNmzNDIkSMlSX/6059UXFwsf38mtAAA3oN3NTdVXl6usLAwc/zdd9+pT58+FiYCAKBlMDPihj777DOziPTs2VO1tbUUEQCA16KMuBHDMHT33XcrOTlZkvTaa6/p+PHjLMsAALwa73JuorS0VB07djTHR48eVWxsrHWBAABoJcyMuIGtW7eaRWTAgAGy2+0UEQCAz6CMWMgwDE2ZMkUpKSmSpDfeeEOHDx9W27ZtLU4GAEDrYZnGImfPnlWXLl3McV5enmJiYixMBACANZgZscDGjRvNIhIfHy+Hw0ERAQD4LMpIKzIMQxMmTNCkSZMkSe+++66ysrLUpg0vAwDAd7FM00pOnTqliIgIc3zixAlFRUVZmAgAAPfAn+St4KOPPjKLyIgRI+R0OikiAAD8C2WkBRmGodGjR+v222+XJK1evVo7d+6Un5+fxckAAHAfLNO0kMLCQkVGRprjoqKiOss0AADgAmZGWsBf/vIXs4jcfPPNcjqdFBEAAOrBzEgzcjqdSkxM1FdffSVJWrdunSZPnmxtKAAA3BxlpJn88MMPdc4Vcvr06TonNQMAAJfGMk0zeOutt8wiMmnSJDmdTooIAAANxMxIEzidTg0ZMkTffPONJGnTpk2aMGGCxakAAPAslJFGOnbsWJ0r6547d8688i4AAGg4lmkaYdmyZWYRueOOO+R0OikiAAA0EjMjLnA4HOrbt6+OHj0qSdq2bZvGjBljcSoAADwbZaSBvv/+e/Xp08ccl5aWKiwszMJEAAB4B5ZpGuCll14yi8h9990nwzAoIgAANBNmRi7DbrcrJiZGhYWFkqT09HQlJSVZnAoAAO9CGanHt99+qwEDBpjjiooKBQcHW5gIAADvxDLNJaSmpppF5MEHH5RhGBQRAABaiMtlJCMjQxMnTlRUVJT8/Py0fv36K+6Tnp6u+Ph4dejQQb169dLy5csbk7XF1dbWKjw8XAsWLJAkffHFF1qyZInFqQAA8G4ul5HKykoNGTJEr7/+eoO2z83N1YQJEzRy5Ejt27dPCxYs0Jw5c7R27VqXw7akgwcPqn379iorK5MkVVVVKTEx0eJUAAB4P5ePGUlJSVFKSkqDt1++fLl69uypl19+WZI0YMAAZWVladGiRZoyZYqrT98i/vjHP+rZZ5+VJD366KNatGiRxYkAAPAdLX4A665duzR27Ng6940bN04rV65UbW2t2rVr97N9bDabbDabOb44W9HcnE6nwsLCVFlZKUnavXu3EhISWuS5AADApbX4AayFhYXq3r17nfu6d+8uu92u06dPX3Kf1NRUhYeHm7eLV8RtbkVFRWYROX/+PEUEAAALtMqnafz8/OqMDcO45P0XzZ8/X6WlpeYtPz+/RXJFRkbK4XDIMAwFBAS0yHMAAIDLa/Flmh49epgnDbuouLhY/v7+6tKlyyX3CQgIaLVy0KYNn24GAMBKLf5OnJiYqG3bttW575NPPlFCQsIljxcBAAC+xeUyUlFRof3792v//v2SLnx0d//+/crLy5N0YYnlnnvuMbd/4IEHdPz4cc2dO1fZ2dl66623tHLlSs2bN695vgMAAODRXF6mycrK0k033WSO586dK0maNm2a3nnnHRUUFJjFRJJiY2O1efNmPfLII1qyZImioqL06quvus3HegEAgLX8jItHk7qxsrIyhYeHq7S0lKvlAgDgIRr6/s3RmwAAwFKUEQAAYCnKCAAAsBRlBAAAWIoyAgAALEUZAQAAlqKMAAAAS1FGAACApSgjAADAUi1+1d7mcPEksWVlZRYnAQAADXXxfftKJ3v3iDJSXl4uSYqJibE4CQAAcFV5ebnCw8Pr/bpHXJvG6XTq5MmTCg0NlZ+fX7M9bllZmWJiYpSfn881b9wcr5Vn4HXyDLxOnsPTXyvDMFReXq6oqCi1aVP/kSEeMTPSpk0bRUdHt9jjh4WFeeSL7It4rTwDr5Nn4HXyHJ78Wl1uRuQiDmAFAACWoowAAABL+XQZCQgI0DPPPKOAgACro+AKeK08A6+TZ+B18hy+8lp5xAGsAADAe/n0zAgAALAeZQQAAFiKMgIAACxFGQEAAJby2TKSkZGhiRMnKioqSn5+flq/fr3VkfATqampGjp0qEJDQxUREaHJkyfryJEjVsfCJSxbtkxxcXHmiZkSExO1ZcsWq2PhClJTU+Xn56eHH37Y6ij4kT/+8Y/y8/Orc+vRo4fVsVqUz5aRyspKDRkyRK+//rrVUVCP9PR0zZw5U5mZmdq2bZvsdrvGjh2ryspKq6PhJ6Kjo/XCCy8oKytLWVlZ+tWvfqVJkybp0KFDVkdDPXbv3q033nhDcXFxVkfBJQwaNEgFBQXm7eDBg1ZHalEecTr4lpCSkqKUlBSrY+Aytm7dWmf89ttvKyIiQnv27FFSUpJFqXApEydOrDN+7rnntGzZMmVmZmrQoEEWpUJ9KioqNHXqVL355ptauHCh1XFwCf7+/l4/G/JjPjszAs9TWloqSercubPFSXA5DodDq1evVmVlpRITE62Og0uYOXOmbrnlFo0ZM8bqKKhHTk6OoqKiFBsbq9/97nc6evSo1ZFalM/OjMCzGIahuXPnasSIEbr22mutjoNLOHjwoBITE3X+/HmFhIRo3bp1GjhwoNWx8BOrV6/W3r17tXv3bqujoB7Dhg3Te++9p759+6qoqEgLFy7UjTfeqEOHDqlLly5Wx2sRlBF4hFmzZunAgQP6/PPPrY6CevTr10/79+9XSUmJ1q5dq2nTpik9PZ1C4kby8/P10EMP6ZNPPlGHDh2sjoN6/PgQgsGDBysxMVHXXHON3n33Xc2dO9fCZC2HMgK3N3v2bG3cuFEZGRmKjo62Og7q0b59e/Xu3VuSlJCQoN27d+uVV17Rn//8Z4uT4aI9e/aouLhY8fHx5n0Oh0MZGRl6/fXXZbPZ1LZtWwsT4lKCg4M1ePBg5eTkWB2lxVBG4LYMw9Ds2bO1bt067dixQ7GxsVZHggsMw5DNZrM6Bn4kOTn5Z5/KuO+++9S/f3898cQTFBE3ZbPZlJ2drZEjR1odpcX4bBmpqKjQ999/b45zc3O1f/9+de7cWT179rQwGS6aOXOm3n//fW3YsEGhoaEqLCyUJIWHhyswMNDidPixBQsWKCUlRTExMSovL9fq1au1Y8eOn30iCtYKDQ392TFXwcHB6tKlC8diuZF58+Zp4sSJ6tmzp4qLi7Vw4UKVlZVp2rRpVkdrMT5bRrKysnTTTTeZ44vrcNOmTdM777xjUSr82LJlyyRJo0ePrnP/22+/rXvvvbf1A6FeRUVFuvvuu1VQUKDw8HDFxcVp69atuvnmm62OBnicH374QXfccYdOnz6tbt26afjw4crMzNRVV11ldbQW42cYhmF1CAAA4Ls4zwgAALAUZQQAAFiKMgIAACxFGQEAAJaijAAAAEtRRgAAgKUoIwAAwFKUEQAAYCnKCAAAsBRlBAAAWIoyAgAALEUZAQAAlvr/ASVpW8xL/64nAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.arrow(1,1,4,3,head_width=0.2)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "ba09306cbb048a98de08b0bf53047f9715d3b35a69d3a22d0c48d222d7ad29fe"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
